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FOREWORD 


This is a supplemental progress report on the research project, “Analysis and Computation 
of Internal Flow-Field in a Scramjet Engine,” for the period ended August 1991. Certain mod- 
ifications were made in the analysis and computational procedures during 1992 after receiving 
inputs from the professional communities at national conferences. The manuscript was prepared, 
in the final form, during the summer of 1993. 

The authors are indebted to Mr. Sutanu Sarkar of the Institute for Computer Applications 
in Science and Engineering (ICASE) at NASA Langley Research Center for his cooperation 
and technical assistance. Partial funding for this research was provided by the NASA Langley 
Research Center through the Grant NAG-1—423. The grant was monitored by Dr. Ajay Kumar 
of Theoretical Flow Physics Branch (Fluid Mechanics Division), Mail Stop 156, NASA Langley 
Research Center, Hampton, Virginia 23681-0001. 



INVESTIGATION OF HIGH-SPEED FREE SHEAR FLOWS USING IMPROVED 
PRESSURE-STRAIN CORRELATED REYNOLDS STRESS TURBULENCE MODEL 

B. Lakshmanan* and S. N. Tiwari 5 
Department of Mechanical Engineering and Mechanics 
Old Dominion University, Norfolk, VA 23529-0247 

ABSTRACT 

A high-speed shear layer is studied using compressibility corrected Reynolds stress turbu- 
lence model which employs newly developed model for pressure-strain correlation. MacCormack 
explicit prediction-corrector method is used for solving the governing equations and the turbu- 
lence transport equations. The stiffeness arising due to source terms in the turbulence equations 
is handled by a semi-implicit numerical technique. Results obtained using the new model show a 
sharper reduction in growth rate with increasing convective Mach number. Some improvements 
were also noted in the prediction of the normalized streamwise stress and Reynolds shear stress. 
The computed results are in good agreement with the experimental data. 
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1. INTRODUCTION 


The recent resurging interest in a High Speed Civil Transport (HSCT) and the National 
Aerospace Plane (NASP) clearly demonstrates the need for advanced propulsion systems for 
supersonic velocities and beyond. Because of the complex nature of the problem, numerous 
research programs have been initiated. One aspect of this research has been directed towards 
detailed understanding of the complex flowfield in the engine over a wide range of operating 
conditions. Computational fluid dynamics (CFD) in conjunction with recent advances in turbu- 
lence modeling is used extensively for detailed investigation of the engine flowfield, as existing 
wind tunnel facilities are inadequate especially in the high Mach number regime. Because of 
the complex nature of the flowfield and increased computational resources necessary, detailed 
simulation of the complete engine problem cannot be considered at the present time. A more 
tractable problem that can be considered in isolation of the complexities introduced by the engine 
geometry is posed by the spatially evolving mixing layer. Currently there is a renewed interest 
in compressible mixing layers for two main reasons. First, mixing layers play an important 
role in many engineering applications, as they are central to many advanced propulsion systems. 
This stems from the fact that they govern the rate of mixing in combustion chambers and are 
also responsible for most of the acoustic noise generated by many propulsion systems. Apart 
from its practical applications, the compressible shear layer problem has remained as the most 
celebrated case for basic testing of transport equation turbulence models. 

It is known that variable density extensions of standard incompresible turbulence models 
[1,2] are inadequate in duplicating the experimentally observed [3,4] reduction in growth rate 
with increasing convective Mach number. This led to attempts by Oh [5], Vandromme [6], and 
Dussage and Quine [7], among others to make modifications to incompressible turbulence models 
in order to obtain successful prediction of the compressible mixing layer. But these modifica- 
tions were developed somewhat on a preliminary fashion and lack the theoretical justifications 
necessary to model the physics of the flow adequately when parameters beyond mean growth 
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rates are considered. Recently, Sarkar et al. [8,9] recognized the importance of compressible 
dissipation and pressure-dilatation which are known to be present in compressible turbulence. 
Simple corrections for compressible dissipation and pressure-dilatation were proposed based on 
direct numerical simulation of compressible isotropic turbulence which can be easily included 
in the existing transport equation turbulence models. In addition, Speziale et al. [10] recently 
developed a new model for pressure-strain correlation which was shown to give improved pre- 
dictions over older models such as the Launder, Reece and Rodi model [2] when applied to 
incompressible, homogeneous turbulent flows. 

During the past decade, considerable progress have been made in the area of computational 
fluid dynamics. Most of the research activities have been centered around computing complex 
three-dimensional flows over realistic aerodynamic configurations. Despite the popularity of the 
existing schemes in computing complex flows, their extension for solving stiff equations are 
limited in the literature, especially within the framework of full Navier-Stokes solvers. 

Recently, Sarkar and Lakshmanan [11] applied the compressibility corrected Reynolds stress 
turbulence model [2] to high speed shear layer using a full Navier-Stokes code. While the 
numerical study was successful, the components of the Reynolds stress were overpredicted. The 
purpose of the study is to incorporate the newly developed pressure-strain correlation in the 
compressibility corrected Reynolds stress model [11]. The resulting turbulence model is then 
applied to the case of high-speed shear layer over a wide range of convective Mach number. 
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2. THEORETICAL FORMULATION 

In this section, essential governing equations, and boundary and initial conditions are 
presented for a compressible supersonic mixing layer evolving downstream of a splitter plate 
(Fig. 2.1). 


2.1 Governing Equations 


The formulation of the problem starts with the Favre-averaged form of the equations 
representing conservation of mass, momentum, energy and turbulence quantities. The overbar is 
used to denote a conventional Reynolds average, whereas the tilde is used to denote the Favre- 
average. For the sake of brevity only the Reynolds stress and dissipation rate transport equations 
are given here. Extensive details of the governing equations and models are given in Ref. [11]. 
Reynolds Stress Equation 


d (p u[u0 + pu k u[uj + T ijk + (<5 ik p"u'' + 5j k p"u|')- 


dt 


( 2 . 1 ) 


(S" k u" + S" k u")] = Pij + Ilij - -p (e s + ec) + - P V' k *j 

V‘dxj Uj dxj 

where P represents production, II is the pressure-strain correlation and T[j k is the diffusive 
transport. The third and fourth term inside the square bracket on the LHS represents the 
contribution due to press ure-velocity and stress-velocity correlations respectively. In the present 
work these correlations are neglected. The third and fourth terms on the RHS represents the 
contribution due to compressible dissipation and pressure-dilatation effects. In the present work, 
these terms are computed using the Sarkar’s model [8,9]. In Eq. (2.1), the production tensor Py 
is defined in an exact manner while models are carried for the remaining terms. 

Production Tensor P;;: 


Pij = -p(u[u' k Uj,k + u'u' k U i)k ) 


( 2 . 2 ) 
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Model for Diffusive Transport: 




Tijk = -C s p 


(u| u j) k + (“X) . + (“!<) 


(2.3) 


where q 12 = u 12 + v' 2 + w' 2 and the model constant C s = 0.018. 
Model for Compressible Dissipation 


l c = ai e 3 Mf 


(2.4) 


Model for Pressure-Dilatation: 

f2 i 


P Hu k,k = 

where M t = 


I q 2 

7 rt 


is the turbulent Mach number and i s is the solenoidal dissipation. Based 


on direct numerical simulation Sarkar [8,9] recommends ai = 0.5, «2 = 0.4 and 03 =0.2. 

Model for Pressure-Strain Correlation 

Two models are employed in the present study for pressure-strain correlation. The first 
model employs the well known Launder, Reece and Rodi model [2] while the second model 
employs the recently developed SSG model [9]. The details of the model are given as follows: 
Original Pressure-Strain Correlation Model [2]: 

n,y = -Cip isbij - c 2 (Pij - ^ s t ^j 

where by is the anisotropy tensor given by 


(2.5) 


bij — 


q 2 


V%} 

3 


(2.6) 


In Eq. (2.6), the model constants are Ci = 3.0 and C 2 = 0.6 
Improved Pressure-Strain Correlation Model [9]: 

1 


*) 


I ^ij — (C3^s *f" C3P )t)jj -f* C 465 

+ (c 5 - ClII 1 ' 2 ) k (s i} - y kk Si^ + Cek ( b ik S jk + b jk s tk (2.7) 


— g ^mnSmn^ij ) ~f* C^k “I - 
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S*j — 2 (“*\i "I" “i.O » J — 2^ u,, - ? 


( 2 . 8 ) 


P — T»y Ujj ( II — bijbij 


(2.9) 


The model constants are 


C 3 = 3.4 , Cl = 1.8 , C 4 = 4.2 

C 5 = 0.8 , C 5 = 1.3 , C 6 = 1.25 , C 7 = 0.4 

The new models for compressible dissipation, pressure-dilatation and pressure-strain correlation 
were chosen because these models have been developed based on direct numerical simulation 
of homogeneous isotropic turbulence. 

Dissipation Transport Equation 

The dissipation transport equation is expressed as 


d d 

a ( ' ,£j + a^ 


pk 7-, &e s 

P u k t s - C ( — u k u i 


r, e s - du i ri P e s 

= -C(l T p U { U-~ C C 2— ~— 

k } dxj k 


( 2 . 10 ) 


The model coefficients in Eq. (2.10) are 


C e i = 1.44 , 0*2 = 1.83, C e = 0.15 


For the present problem, we need to solve the Navier-Stokes equations along with the 
equations of state, to obtain the mean variables p, u, v and E. In the case of the plane shear 
layer, the Reynolds stress tensor has four nonzero component: u' 2 , u ,2 ) w' 2 and u'v', which are 
solved by the corresponding components of Eq. (2.1). The equation for the solinoidal dissipation 
e s completes the set of governing equations. Thus, a system of nine coupled, nonlinear, partial 
differential equations along with an appropriate set of initial and boundary conditions must be 


solved. 
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2.2 Boundary and Initial Conditions 

Since the governing equations are elliptic in nature, the boundary conditions have to be 
specified along all four boundaries. These include inflow, outflow, and outer boundaries Gower 
and upper boundaries), respectively. At the inflow boundary (at = 0.0), profiles are specified for 
the velocities, static pressure, static temperature, turbulent stresses, and the turbulent dissipation 
rate. Since we are interested in the downstream fully develop regime, the specific form of the 
inlet profiles is not crucial. 

The outer boundaries always remain in the ffeestream, and the appropriate boundary condi- 
tion is to assume that the normal derivative of the flow variables vanish along those boundaries. 
The gradient boundary conditions not only preserve the freestream values along the outer bound- 
aries but also provide nonreflective conditions for the outgoing waves. The outflow boundary 
(jc = Xmax) is always supersonic and, hence, the values of mean flow and turbulence quantities 
are determined by zeroth-order extrapolation from upstream values. Along with the boundary 
conditions, the governing equations also require a set of initial conditions. The initial conditions 
at time t = 0 for all of the variables are obtained by simply propagating the inflow profiles 
throughout the computational domain. Having specified all of the boundary and initial data, 
the equations are marched in time until the residual based on pu decreases by six orders of 
magnitude, indicating that a converged solution has been obtained. 
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3. METHOD OF SOLUTION 


The transport equations for the mean flow and Reynolds stresses are written in the physical 
domain and must be transformed to the computational domain using an appropriate coordinate 
transformation. For the physical problem under consideration, an algebraic grid generation 
technique is used to generate the mesh. In the physical domain, a uniform grid is used in the 
axial direction and in the normal direction the grid lines are clustered near regions where strong 
gradients exist. A uniform mesh is used in the computational domain. The governing equations 
are first cast into a vector form, where U is the dependent variable vector consisting of nine 
components; the vectors F and G, respectively, denote the x and y destruction and redistribution 
of the Reynolds stresses. To numerically obtain the solution for the vector U, the governing 
equations are then transformed from the physical domain to the computational domain, giving 
the following system of equations, 


dU_ dF_ dG_£r 
dt + d£ + dy ~ 


(3.1) 


where 


0 = JU, H = JH 

F = Fy v - GX Vl G = Gx^ - Fy ^ J = x^y v - y^x v 

In Eq. (3.1), a superscript Q denotes quantities in the transformed system, (a:^, x v , y^,y v ) represent 
the metrics of the transformation, and / denotes the Jacobian of the transformation. If the physical 
grid is given, the metrics and the Jacobian of the transformation can be computed easily. 

The governing equations are integrated explicitly in time using the unsplit MacCormack 
predictor-corrector scheme [12]. During a specific numerical sweep, the inviscid fluxes and 
the first-derivative terms in the source vector H are backward differenced in the predictor step 
and forward differenced in the corrector step. Second-order central differences are used for the 
viscous and heat flux terms. Hence, the complete scheme for both the predictor and corrector 
steps can be expressed as follows 
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Predictor: 


V„G", „ 

At/"/ 1 = -At -4^- + -4-^ - # r, 


(3.2a) 


v°r = i>Zi + 


(3.2b) 


Corrector: 


At/"/ 1 = - At A ffi + A ^ ,,j - ^" +i 

l A£ A TJ l ' } 


(3.3a) 


J7// 1 = 1/2 (f// + L// 1 + At/// 1 ) 


(3.3b) 


The composite numerical scheme is second-order accurate in both time and space and, being 
an explicit scheme, is conditionally restricted by the Courant and viscous stability limits of the 
governing equations. The solution procedure requires no scalar or block tridiagonal inversions. 
The flow- field is advanced from time level n to n +1 and this process is continued until the 
desired integration time or steady state has been reached. Since the Reynolds stress transport 
equations contain stiff source terms, the maximum Courant-Friedricks-Lewy (CFL) number used 
in the computation was limited to 0.5. 

The numerical code used in this study is a two-dimensional, Navier-Stokes solver [13] written 
in a generalize body-oriented coordinate system. As such, various two-dimensional free shear 
flows and wall bounded flows can be handled by the numerical code. The code in its original 
form used a second-order spatially and temporally accurate two-step MacCormack scheme. The 
later versions of the code employ a variety of higher order compact algorithms [14] (fourth and 
sixth order) and various upwind schemes. Local time stepping and residual smoothing options 
are also available in the code to accelerate the convergence to steady state. In the present research 
work, the capabilities of the code are further enhanced by adding a second-order Reynolds stress 
model as a turbulence closure. 
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4. RESULTS AND DISCUSSION 

A shear layer developing due to mixing between parallel supersonic streams in chosen as a 
test case to validate the newly developed turbulence model. Extensive results have been obtained 
over a wide range of convective Mach number. These results serve two purposes. First, the 
numerical procedure and the developed computer code are validated by comparing with the 
experimental results [4, 15-18]. Secondly, they explain the characteristic behavior of high-speed 
shear layers. It should be pointed out that all the results reported in this paper are computed 
without artificial dissipation added to the numerical scheme. 

Results were obtained for extensive range of convective Mach number by varying the upper 
speed Ui from 900 m/sec to 4000 m/sec while maintaining the lower speed U 2 at 800 m/sec. 
The static conditions for the two incident streams were assumed to be equal. Specifically the 
results were obtained for the following conditions: 

Poo “ 101325 N/m 2 , = 800 K, Poo = 0.44 kg/m 3 


Computations were carried out for a shear layer (Fig. 2.1) developing over a length of 10 
cm. The height of the domain was assumed to be 5 cm. Before discussing the results, a few 
definitions are in order. It is well known from the experiments that a fully-developed shear layer 
spreads linearly, and that the growth rate dSIdx can be expressed by the following relation: 


dS_ _ / ui - U2 

dx \lil + U 2 


(4.1) 


where 6(x) denotes the thickness of the shear layer, and is approximately constant. The 
shear layer thickness <5(x) represents the distance between the two cross-stream locations where 
the normalized streamwise velocity u* = (u — u 2 )/(u\ — U 2 ) is 0.1 and 0.9 respectively. The 
convective Mach number M c is defined by the following relationship: 


M c = 


/til — U2 
\ a l + a 2 


(4.2) 


where aj and a 2 are the speeds of sound in the two layers. 
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The calculations were carried out using 201x51 mesh. A uniformly spaced grid was used 
in the streamwise direction while a stretched grid was employed in the cross-stream direction. 
A grid refinement study was carried out in the cross-stream direction for all convective Mach 
numbers using 101 and 201 points. Grid independent results were obtained between 101 and 201 
points and the variation in the solution between 51 and 101 points was found to be less than 1 %. 

Figures 4.1 - 4.4 show the results obtained using the compressibility corrected form of the 
original model and the improved model for a particular set of conditions of the shear layer 
having primary stream velocity u\ = 2500 m / s - 

The velocity and the Reynolds stress profiles were plotted as a function of the similarity 
variable 77 = (y — y c )/8 where y is the normal coordinate and y c is the normal coordinate 
location where u*=0.5. 

Figure 4.1 shows that the mixing layer thickness 8(x) increases linearly after an initial 
development phase. It is evident that at the outflow boundary of the computational domain, the 
linearly growing regime is well established and the mean velocity profiles have achieved self 
similar form. The improved model predicts a slightly reduced mixing layer thickness compared 
to the predictions shown by the original model. 

Figures 4.2 - 4.4 show the comparison of the fully developed mean velocity and turbulent 
stress profiles computed by the two models. These profiles clearly display the asymmetry present 
in the flow field as shown by the downward movement of the center of the mixing layer (y 
coordinate location where u*=0.5) with axial location (Fig. 4.5). The improved model yields 
slightly reduced peak values for the normal stress and Reynolds shear Stress. These results clearly 
indicate a greater-penetration of the flow into the low-speed than the corresponding penetration 
into the upper high-speed side of the domain. 

Figure 4.6 shows the normalized growth rate for a fully developed flow ((C ^) 0 being the 
incompressible value which was obtained by calculating case with a small M c ) as a function of 
convective Mach number. It is seen that without compressibility effects (a=0.0) both the models 
show only a mild decrease in growth rate with increasing convective Mach number. However, 
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inclusion of the compressibility correction to the models ( 0 = 1 . 0 ) lead to a sharper reduction in 
growth rate with increasing convective Mach number. Improved results for the growth rate are 
obtained using the compressibility corrected form of the new model as compared to the original 
model especially in the high convective Mach number regime. Figures 4.7 and 4.8 show the 
computed peak normal stress and Reynolds shear Stress and comparison with the experimental 
data [4, 16-18]. The peak values of the normalized stress components decrease significantly 
with increasing convective Mach number. This is thought to be a direct consequence of the 
compressibility correction which has a dissipative effect on the growth rate of the turbulent 
kinetic energy and shear stress. It is clear that results obtained using the improved model are 
well within the range of the experimental data unlike those of the original model. 
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5. CONCLUSIONS 

Initially, a second-order turbulence closure (employing Lauder, Reece and Rodi model for 
pressure-strain correlation) without any compressibility correction was applied to high-speed 
shear layer. The results confirmed that variable density extensions of incompressible turbulence 
models were inadequate in duplicating the experimentally reduction in growth rate with increasing 
convective Mach number. When compressibility effects were included in the models the results 
showed a dramatic reduction in growth rate and peak values of the normalized Reynolds stress 
components with increasing convective Mach number. In addition, the results obtained using the 
new model for pressure-strain correlation gave improved agreement with the experimental data 
for normalized growth rate and peak values of the Reynolds stress components with increasing 
convective Mach number. 
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Fig. 4.1 Axial variation of Hip shear layer thickness 
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Fig. 4.2 Comparison of axial velocity profile at x = 10 cm 
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Fig. 4.4 Comparison of I hr Reynolds shear stress profile at x = 10 cm 



24.2 



x, cm 


Fig. 4.5 Variation of the y coordinate location where u = 0.5 with axial 
location 
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Fig. 4.7 Variation of Hip maximum streamwlse stress with the convective 
Mach number 
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Fig. 4.8 Variation of tlir maximum Reynold shear stress with the 
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